Modelo híbrido ARIMA - XGBoost para predecir el Precio del Gas Natural#
ARIMA (Autoregressive Integrated Moving Average), es un modelo estadístico clásico utilizado para analizar y predecir series de tiempo estacionarias o que se pueden hacer estacionarias mediante diferenciación. ARIMA(p,d,q) tiene tres componentes principales:
Autorregresivo (AR): Utiliza valores pasados de la serie para predecir valores futuros. p es el orden del modelo AR, que representa el número de valores pasados utilizados.
Integrado (I): La diferenciación (d) se aplica a la serie para lograr la estacionariedad, es decir, eliminar tendencias y estacionalidad.
Media móvil (MA): Utiliza los errores pasados (residuos) de la serie para predecir valores futuros. q es el orden del modelo MA, que indica el número de errores pasados.
ARIMA es eficaz para predecir series de tiempo con tendencias y patrones estacionales claros. Su fortaleza reside en su base teórica sólida y su capacidad para capturar relaciones lineales entre valores pasados y futuros. Sin embargo, ARIMA puede tener limitaciones al modelar no linealidades complejas.
Carga y Procesamiento de los datos#
Antes de desarrollar los modelos, comenzamos importando las librerías requeridas y cargando la base de datos de precios de cierre del gas natural. Posteriormente convertimos el formato de la columna ‘date’ a fecha y hora, la establecemos como índice, convertimos las columnas a numerica y se eliminan las filas que contienen valores NaN, lo cual es necesario porque el desplazamiento crea NaN en las primeras filas.
import pandas as pd
import numpy as np
from statsmodels.tsa.arima.model import ARIMA
from sklearn.model_selection import TimeSeriesSplit, RandomizedSearchCV, train_test_split
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error
import matplotlib.pyplot as plt
import plotly.io as pio
import plotly.express as px
import plotly.offline as py
pio.renderers.default = "notebook"
import warnings
warnings.filterwarnings("ignore")
data=pd.read_csv(r"C:\Users\claud\OneDrive - Universidad del Norte\Escritorio\Series de tiempo\Natural_Gas_data.csv")
data['date'] = pd.to_datetime(data['date'])
data.set_index('date', inplace=True)
data_index = data.index
Partición de los datos en entrenamiento, prueba y test
y = data['close'].values
# Primero, dividimos en train + val y test (85% para train + val, 15% para test)
y_train_val, y_test = train_test_split(y, test_size=0.15, shuffle=False)
# Luego, dividimos el train + val en train (70%) y val (30% de ese subconjunto)
y_train, y_val = train_test_split(y_train_val, test_size=0.3, shuffle=False)
# Verificando las formas de los conjuntos
print(f'Tamaño de entrenamiento: {len(y_train)}')
print(f'Tamaño de validación: {len(y_val)}')
print(f'Tamaño de prueba: {len(y_test)}')
Tamaño de entrenamiento: 3558
Tamaño de validación: 1525
Tamaño de prueba: 897
data_index_full = data_index[:len(y_train) + len(y_val) + len(y_test)]
# Extender las listas para alinear los datos con el índice completo
y_train_extended = np.concatenate([y_train, [None] * (len(y_val) + len(y_test))])
y_val_extended = np.concatenate([[None] * len(y_train), y_val, [None] * len(y_test)])
y_test_extended = np.concatenate([[None] * (len(y_train) + len(y_val)), y_test])
# Crear DataFrame para la gráfica
df_plot = pd.DataFrame({
'Fecha': data_index_full,
'Datos de entrenamiento': y_train_extended,
'Datos de validación': y_val_extended,
'Datos de prueba': y_test_extended,
})
# Graficar con Plotly
fig = px.line(
df_plot, x='Fecha',
y=[
'Datos de entrenamiento',
'Datos de validación',
'Datos de prueba'
],
title="Particiones serie: precio de cierre gas natural", labels={'value': 'Valor', 'Fecha': 'Fecha'}
)
fig.show()
def calculate_and_store_metrics(model, model_name, y_train, y_train_pred, y_val, y_val_pred, y_test, y_pred):
# Cálculo de métricas para el conjunto de prueba
mape_test = mean_absolute_percentage_error(y_test, y_pred)
rmse_test = mean_squared_error(y_test, y_pred, squared=False)
# Cálculo de métricas para el conjunto de validación
mape_val = mean_absolute_percentage_error(y_val, y_val_pred)
rmse_val = mean_squared_error(y_val, y_val_pred, squared=False)
# Almacenando las métricas en la lista
metrics_list.append({
'Modelo': model_name,
'MAPE Prueba': mape_test,
'RMSE Prueba': rmse_test,
'MAPE Validación': mape_val,
'RMSE Validación': rmse_val
})
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy.stats import shapiro # Corrección: Usamos shapiro desde scipy
def check_model_assumptions(model, y_train, y_train_pred, y_val, y_val_pred, y_test, y_test_pred):
# Calcular los residuales
resid_train = y_train - y_train_pred
resid_val = y_val - y_val_pred
resid_test = y_test - y_test_pred
print("\n----- Pruebas de Hipótesis sobre los Supuestos -----")
# Prueba de Normalidad (Shapiro-Wilk) para entrenamiento
stat, p_value = shapiro(resid_train)
print(f"Shapiro-Wilk p-value (Entrenamiento): {p_value:.4f}")
print("Conclusión:", "NO sigue una distribución normal.\n" if p_value < 0.05 else "Sigue una distribución normal.\n")
# Prueba de Autocorrelación (Ljung-Box) para entrenamiento
lb_test = sm.stats.acorr_ljungbox(resid_train, lags=[10], return_df=True)
lb_p_value = lb_test['lb_pvalue'].values[0]
print(f"Ljung-Box p-value (Entrenamiento): {lb_p_value:.4f}")
print("Conclusión:", "Autocorrelación significativa.\n" if lb_p_value < 0.05 else "Sin autocorrelación significativa.\n")
# Prueba de Homocedasticidad (Breusch-Pagan) para entrenamiento
exog = sm.add_constant(np.arange(len(resid_train)).reshape(-1, 1))
_, bp_p_value, _, _ = sm.stats.het_breuschpagan(resid_train, exog)
print(f"Breusch-Pagan p-value (Entrenamiento): {bp_p_value:.4f}")
print("Conclusión:", "Varianza no constante.\n" if bp_p_value < 0.05 else "Varianza constante.\n")
# Crear figuras en Matplotlib
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
fig.suptitle('Verificación de los Supuestos del Modelo')
# Histograma de los residuales de entrenamiento
axes[0, 0].hist(resid_train, bins=20, edgecolor='black')
axes[0, 0].set_title('Histograma (Entrenamiento)')
axes[0, 0].set_xlabel('Residual')
axes[0, 0].set_ylabel('Frecuencia')
# Q-Q Plot de los residuales de entrenamiento
sm.qqplot(resid_train, line='s', ax=axes[0, 1])
axes[0, 1].set_title('Q-Q Plot (Entrenamiento)')
# ACF de los residuales de entrenamiento
sm.graphics.tsa.plot_acf(resid_train, ax=axes[0, 2], lags=30)
axes[0, 2].set_title('ACF (Entrenamiento)')
# Histograma de los residuales de validación
axes[1, 0].hist(resid_val, bins=20, edgecolor='black')
axes[1, 0].set_title('Histograma (Validación)')
axes[1, 0].set_xlabel('Residual')
axes[1, 0].set_ylabel('Frecuencia')
# Q-Q Plot de los residuales de validación
sm.qqplot(resid_val, line='s', ax=axes[1, 1])
axes[1, 1].set_title('Q-Q Plot (Validación)')
# ACF de los residuales de validación
sm.graphics.tsa.plot_acf(resid_val, ax=axes[1, 2], lags=30)
axes[1, 2].set_title('ACF (Validación)')
plt.tight_layout(rect=[0, 0, 1, 0.96]) # Ajustar layout para evitar superposición
plt.show()
def plot_model_results_plotly(y_train, y_train_pred, y_val, y_val_pred, y_test, y_pred, data_index, train_size, val_size, title):
# Crear un índice completo basado en la longitud total de los datos
data_index_full = data_index[:len(y_train) + len(y_val) + len(y_test)]
# Extender las listas para alinear los datos con el índice completo
y_train_extended = np.concatenate([y_train, [None] * (len(y_val) + len(y_test))])
y_train_pred_full = np.concatenate([y_train_pred, [None] * (len(y_train_extended) - len(y_train_pred))])
y_val_extended = np.concatenate([[None] * len(y_train), y_val, [None] * len(y_test)])
y_val_pred_full = np.concatenate([[None] * len(y_train), y_val_pred, [None] * (len(y_test))])
y_test_extended = np.concatenate([[None] * (len(y_train) + len(y_val)), y_test])
y_pred_full = np.concatenate([[None] * (len(y_train) + len(y_val)), y_pred])
# Crear DataFrame para la gráfica
df_plot = pd.DataFrame({
'Fecha': data_index_full,
'Datos de entrenamiento': y_train_extended,
'Modelo entrenado': y_train_pred_full,
'Datos de validación': y_val_extended,
'Predicción validación': y_val_pred_full,
'Datos de prueba': y_test_extended,
'Predicción prueba': y_pred_full
})
# Graficar con Plotly
fig = px.line(
df_plot, x='Fecha',
y=[
'Datos de entrenamiento', 'Modelo entrenado',
'Datos de validación', 'Predicción validación',
'Datos de prueba', 'Predicción prueba'
],
title=title, labels={'value': 'Valor', 'Fecha': 'Fecha'}
)
fig.show()
Modelo ARIMA#
Luego, se crea un modelo ARIMA con el orden (5, 1, 0) usando la serie de precios de cierre. El orden (p, d, q) representa el orden autorregresivo (p), el grado de diferenciación (d) y el orden de media móvil (q). (5,1,0) indica un modelo AR(5) con una diferenciación de primer orden.
El modelo se ajusta a los datos usando model_arima.fit(), se generan predicciones para toda la serie temporal y se calculan los residuos restando las predicciones del modelo ARIMA de los valores reales.
def rolling_forecast(order, y_train, y_val, y_test):
history = list(y_train) # Usar datos de entrenamiento como historial inicial
val_predictions = []
test_predictions = []
# Predecir para el conjunto de validación
for t in range(len(y_val)):
model = ARIMA(history, order=order).fit()
yhat = model.forecast()[0]
val_predictions.append(yhat)
history.append(y_val[t]) # Actualizar el historial con el dato real
# Predecir para el conjunto de prueba
for t in range(len(y_test)):
model = ARIMA(history, order=order).fit()
yhat = model.forecast()[0]
test_predictions.append(yhat)
history.append(y_test[t]) # Actualizar el historial con el dato real
return np.array(val_predictions), np.array(test_predictions)
def optimize_model(order_range, y_train):
"""Optimiza el modelo ARIMA para los parámetros dados."""
best_aic = float("inf")
best_params = None
best_model = None
for p in order_range[0]:
for d in order_range[1]:
for q in order_range[2]:
try:
model = ARIMA(y_train, order=(p, d, q)).fit()
aic = model.aic
if aic < best_aic:
best_aic = aic
best_params = (p, d, q)
best_model = model
except Exception as e:
print(f"Error: {e}")
continue
print(f"Mejores parámetros: p={best_params[0]}, d={best_params[1]}, q={best_params[2]} con AIC={best_aic}")
return best_model, best_params
def evaluate_models(y_train, y_val, y_test, data):
"""Evalúa todos los modelos y selecciona el mejor."""
models = [
('AR', (p_values, [0], [0])),
('MA', ([0], [0], q_values)),
('ARIMA', (p_values, d_values, q_values))
]
best_aic = float("inf")
best_model_name = None
best_model = None
best_params = None
for model_name, order_range in models:
print(f"Evaluando: {model_name}")
model, params = optimize_model(order_range, y_train)
if model.aic < best_aic:
best_aic = model.aic
best_model_name = model_name
best_model = model
best_params = params
print(f"\nEl mejor modelo es: {best_model_name} con AIC={best_aic}")
# Rolling forecast para validación y prueba
y_val_pred, y_test_pred = rolling_forecast(best_params, y_train, y_val, y_test)
y_train_pred = best_model.fittedvalues
# Calcular y almacenar métricas
calculate_and_store_metrics(best_model, best_model_name, y_train, y_train_pred, y_val, y_val_pred, y_test, y_test_pred)
# Graficar los resultados
plot_model_results_plotly(y_train, y_train_pred, y_val, y_val_pred, y_test, y_test_pred, data.index, len(y_train), len(y_val), best_model_name)
check_model_assumptions(best_model, y_train, y_train_pred, y_val, y_val_pred, y_test, y_test_pred)
return(best_model, best_model_name, y_train_pred, y_val_pred, y_test_pred)
metrics_list=[]
p_values = range(0, 3)
d_values = range(0, 2)
q_values = range(0, 6)
best_model, best_model_name, y_train_pred, y_val_pred, y_test_pred = evaluate_models(y_train, y_val, y_test, data)
Evaluando: AR
Mejores parámetros: p=1, d=0, q=0 con AIC=-911.4617908954825
Evaluando: MA
Mejores parámetros: p=0, d=0, q=5 con AIC=3438.3500291627465
Evaluando: ARIMA
Mejores parámetros: p=2, d=1, q=5 con AIC=-923.2849978518195
El mejor modelo es: ARIMA con AIC=-923.2849978518195
----- Pruebas de Hipótesis sobre los Supuestos -----
Shapiro-Wilk p-value (Entrenamiento): 0.0000
Conclusión: NO sigue una distribución normal.
Ljung-Box p-value (Entrenamiento): 0.9600
Conclusión: Sin autocorrelación significativa.
Breusch-Pagan p-value (Entrenamiento): 0.0000
Conclusión: Varianza no constante.
Ahora tomamos los residuos del modelo para hacer el modelo de redes neuronales
resid_train = y_train - y_train_pred
resid_val = y_val - y_val_pred
resid_test = y_test - y_test_pred
#g residuos
data_index_full = data_index[:len(resid_train) + len(resid_val) + len(resid_test)]
# Extender las listas para alinear los datos con el índice completo
y_train_extended = np.concatenate([resid_train, [None] * (len(resid_val) + len(resid_test))])
y_val_extended = np.concatenate([[None] * len(resid_train), resid_val, [None] * len(resid_test)])
y_test_extended = np.concatenate([[None] * (len(resid_train) + len(resid_val)), resid_test])
# Crear DataFrame para la gráfica
df_plot = pd.DataFrame({
'Fecha': data_index_full,
'resid de entrenamiento': y_train_extended,
'resid de validación': y_val_extended,
'resid de prueba': y_test_extended,
})
# Graficar con Plotly
fig = px.line(
df_plot, x='Fecha',
y=[
'resid de entrenamiento',
'resid de validación',
'resid de prueba'
],
title="Particiones serie: precio de cierre gas natural", labels={'value': 'Valor', 'Fecha': 'Fecha'}
)
fig.show()
import numpy as np
import pandas as pd
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, Dropout
from tensorflow.keras.callbacks import EarlyStopping
from sklearn.metrics import mean_squared_error
# Paso 1: Función para crear X e y usando ventanas deslizantes
def makeXy(ts, nb_timesteps):
"""
Transforma una serie de tiempo en conjuntos X (regresores) e y (objetivo) usando ventanas deslizantes.
"""
X = []
y = []
for i in range(nb_timesteps, ts.shape[0]):
X.append(list(ts[i - nb_timesteps:i])) # Ventana deslizante como regresores
y.append(ts[i]) # El siguiente valor como objetivo
X, y = np.array(X), np.array(y)
return X, y
# Paso 2: Función para crear el modelo MLP
def create_mlp_model(input_shape, neurons, dropout_rate):
model = Sequential()
model.add(Dense(units=neurons[0], activation='relu', input_shape=(input_shape,)))
model.add(Dropout(dropout_rate))
model.add(Dense(units=neurons[1], activation='relu'))
model.add(Dropout(dropout_rate))
model.add(Dense(units=1)) # Capa de salida
model.compile(optimizer='adam', loss='mean_squared_error')
return model
# Paso 3: Función para entrenar y evaluar modelos con diferentes combinaciones de hiperparámetros
def train_and_evaluate_models(X_train, y_train, X_val, y_val, neuron_combinations, dropout_combinations, batch_sizes, epochs=50):
results = []
best_model = None
best_val_loss = float('inf')
best_params = {}
for neurons in neuron_combinations:
for dropout_rate in dropout_combinations:
for batch_size in batch_sizes:
# Crear el modelo con la configuración actual
model = create_mlp_model(X_train.shape[1], neurons, dropout_rate)
# Entrenar el modelo con EarlyStopping
history = model.fit(
X_train, y_train,
epochs=epochs,
batch_size=batch_size,
validation_data=(X_val, y_val),
verbose=0,
callbacks=[EarlyStopping(patience=10, restore_best_weights=True)]
)
# Evaluar el modelo en los datos de validación
val_loss = model.evaluate(X_val, y_val, verbose=0)
# Guardar los resultados de esta combinación
results.append({
'neurons': neurons,
'dropout_rate': dropout_rate,
'batch_size': batch_size,
'val_loss': val_loss
})
# Verificar si este modelo tiene la mejor pérdida de validación
if val_loss < best_val_loss:
best_val_loss = val_loss
best_model = model
best_params = {
'neurons': neurons,
'dropout_rate': dropout_rate,
'batch_size': batch_size
}
# Imprimir el progreso
print(f"Neurons: {neurons}, Dropout: {dropout_rate}, Batch Size: {batch_size}, Validation Loss: {val_loss}")
# Convertir resultados a un DataFrame para una mejor visualización
results_df = pd.DataFrame(results)
print("\nMejor modelo:")
print(f"Neurons: {best_params['neurons']}, Dropout: {best_params['dropout_rate']}, Batch Size: {best_params['batch_size']}, Validation Loss: {best_val_loss}")
return best_model, results_df, best_params
# Ejemplo de uso de la función
# Supongamos que resid_train, resid_val, resid_test están definidos como series de tiempo de numpy
# Configuración de parámetros
nb_timesteps = 5 # Número de pasos de tiempo en cada ventana deslizante
neuron_combinations = [(64, 32), (128, 64), (32, 16)]
dropout_combinations = [0.2, 0.3, 0.4]
batch_sizes = [16, 32, 64]
# Crear X e y para train, val y test usando la función makeXy
X_train_resd, y_train_resd = makeXy(resid_train, nb_timesteps)
X_val_resd, y_val_resd = makeXy(resid_val, nb_timesteps)
X_test_resd, y_test_resd = makeXy(resid_test, nb_timesteps)
# Entrenar el modelo y obtener los mejores parámetros
best_model_resd, results_df_resd, best_params_resd = train_and_evaluate_models(X_train_resd, y_train_resd, X_val_resd, y_val_resd, neuron_combinations, dropout_combinations, batch_sizes)
Neurons: (64, 32), Dropout: 0.2, Batch Size: 16, Validation Loss: 0.007824108935892582
Neurons: (64, 32), Dropout: 0.2, Batch Size: 32, Validation Loss: 0.007852353155612946
Neurons: (64, 32), Dropout: 0.2, Batch Size: 64, Validation Loss: 0.00779171846807003
Neurons: (64, 32), Dropout: 0.3, Batch Size: 16, Validation Loss: 0.007849524728953838
Neurons: (64, 32), Dropout: 0.3, Batch Size: 32, Validation Loss: 0.00781164038926363
Neurons: (64, 32), Dropout: 0.3, Batch Size: 64, Validation Loss: 0.007905117236077785
Neurons: (64, 32), Dropout: 0.4, Batch Size: 16, Validation Loss: 0.007881257683038712
Neurons: (64, 32), Dropout: 0.4, Batch Size: 32, Validation Loss: 0.007807622663676739
Neurons: (64, 32), Dropout: 0.4, Batch Size: 64, Validation Loss: 0.00783879030495882
Neurons: (128, 64), Dropout: 0.2, Batch Size: 16, Validation Loss: 0.007814438082277775
Neurons: (128, 64), Dropout: 0.2, Batch Size: 32, Validation Loss: 0.007789664436131716
Neurons: (128, 64), Dropout: 0.2, Batch Size: 64, Validation Loss: 0.007861248217523098
Neurons: (128, 64), Dropout: 0.3, Batch Size: 16, Validation Loss: 0.007826651446521282
Neurons: (128, 64), Dropout: 0.3, Batch Size: 32, Validation Loss: 0.00773252546787262
Neurons: (128, 64), Dropout: 0.3, Batch Size: 64, Validation Loss: 0.007812670432031155
Neurons: (128, 64), Dropout: 0.4, Batch Size: 16, Validation Loss: 0.007832692004740238
Neurons: (128, 64), Dropout: 0.4, Batch Size: 32, Validation Loss: 0.007816466502845287
Neurons: (128, 64), Dropout: 0.4, Batch Size: 64, Validation Loss: 0.007841799408197403
Neurons: (32, 16), Dropout: 0.2, Batch Size: 16, Validation Loss: 0.007831436581909657
Neurons: (32, 16), Dropout: 0.2, Batch Size: 32, Validation Loss: 0.00782731268554926
Neurons: (32, 16), Dropout: 0.2, Batch Size: 64, Validation Loss: 0.007894482463598251
Neurons: (32, 16), Dropout: 0.3, Batch Size: 16, Validation Loss: 0.007827219553291798
Neurons: (32, 16), Dropout: 0.3, Batch Size: 32, Validation Loss: 0.0077741555869579315
Neurons: (32, 16), Dropout: 0.3, Batch Size: 64, Validation Loss: 0.007820922881364822
Neurons: (32, 16), Dropout: 0.4, Batch Size: 16, Validation Loss: 0.007872060872614384
Neurons: (32, 16), Dropout: 0.4, Batch Size: 32, Validation Loss: 0.00786652509123087
Neurons: (32, 16), Dropout: 0.4, Batch Size: 64, Validation Loss: 0.007905959151685238
Mejor modelo:
Neurons: (128, 64), Dropout: 0.3, Batch Size: 32, Validation Loss: 0.00773252546787262
best_model_resd.summary()
Model: "sequential_67"
_________________________________________________________________
Layer (type) Output Shape Param #
=================================================================
dense_201 (Dense) (None, 128) 768
dropout_134 (Dropout) (None, 128) 0
dense_202 (Dense) (None, 64) 8256
dropout_135 (Dropout) (None, 64) 0
dense_203 (Dense) (None, 1) 65
=================================================================
Total params: 9,089
Trainable params: 9,089
Non-trainable params: 0
_________________________________________________________________
# Predicciones en train, val y test
train_predictions_resd = best_model_resd.predict(X_train_resd).ravel()
val_predictions_resd = best_model_resd.predict(X_val_resd).ravel()
test_predictions_resd = best_model_resd.predict(X_test_resd).ravel()
112/112 [==============================] - 0s 1ms/step
48/48 [==============================] - 0s 822us/step
28/28 [==============================] - 0s 880us/step
# Calcular el MSE en cada conjunto
train_mse_resd = mean_squared_error(y_train_resd, train_predictions_resd)
val_mse_resd = mean_squared_error(y_val_resd, val_predictions_resd)
test_mse_resd = mean_squared_error(y_test_resd, test_predictions_resd)
print("\nMSE en conjunto de entrenamiento:", train_mse_resd)
print("MSE en conjunto de validación:", val_mse_resd)
print("MSE en conjunto de prueba:", test_mse_resd)
MSE en conjunto de entrenamiento: 0.04327132513049912
MSE en conjunto de validación: 0.007732527109247361
MSE en conjunto de prueba: 0.05698560495812593
# Mostrar los resultados ordenados por la pérdida de validación
results_df = results_df_resd.sort_values(by='val_loss', ascending=True)
print(results_df)
neurons dropout_rate batch_size val_loss
13 (128, 64) 0.3 32 0.007733
22 (32, 16) 0.3 32 0.007774
10 (128, 64) 0.2 32 0.007790
2 (64, 32) 0.2 64 0.007792
7 (64, 32) 0.4 32 0.007808
4 (64, 32) 0.3 32 0.007812
14 (128, 64) 0.3 64 0.007813
9 (128, 64) 0.2 16 0.007814
16 (128, 64) 0.4 32 0.007816
23 (32, 16) 0.3 64 0.007821
0 (64, 32) 0.2 16 0.007824
12 (128, 64) 0.3 16 0.007827
21 (32, 16) 0.3 16 0.007827
19 (32, 16) 0.2 32 0.007827
18 (32, 16) 0.2 16 0.007831
15 (128, 64) 0.4 16 0.007833
8 (64, 32) 0.4 64 0.007839
17 (128, 64) 0.4 64 0.007842
3 (64, 32) 0.3 16 0.007850
1 (64, 32) 0.2 32 0.007852
11 (128, 64) 0.2 64 0.007861
25 (32, 16) 0.4 32 0.007867
24 (32, 16) 0.4 16 0.007872
6 (64, 32) 0.4 16 0.007881
20 (32, 16) 0.2 64 0.007894
5 (64, 32) 0.3 64 0.007905
26 (32, 16) 0.4 64 0.007906
# Predicciones en train, val y test
plot_model_results_plotly(y_train_resd, train_predictions_resd, y_val_resd, val_predictions_resd, y_test_resd, test_predictions_resd, data.index, len(y_train_resd), len(y_val_resd), 'Modelo híbrido')
Sumamos los residuos prredichos a las predicciones arima
a las prediicones de arima se le quitan los primero 5 datos para ajustarlo con las predicciones de los residuos
train_join=train_predictions_resd + y_train_pred[5:]
val_join=val_predictions_resd + y_val_pred[5:]
test_join= test_predictions_resd+ y_test_pred[5:]
plot_model_results_plotly(y_train[5:],train_join , y_val[5:], val_join, y_test[5:], test_join, data.index, len(y_train), len(y_val), 'predicciones parciales')
Calcualmos los residuos del modelo híbrido
resid_join_train= y_train[5:]- train_join
resid_join_val= y_val[5:]- val_join
resid_join_test= y_test[5:]-test_join
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error
from scipy.stats import shapiro, ttest_1samp
from statsmodels.stats.diagnostic import acorr_ljungbox, het_breuschpagan
import statsmodels.api as sm
# Función para analizar residuales
def analyze_residuals(residuals, y_true, predictions, dataset_name=""):
# Calcular MAPE y MSE
mape = mean_absolute_percentage_error(y_true[5:], predictions) # Ignorando los primeros 5 valores como en los residuales
mse = mean_squared_error(y_true[5:], predictions)
# Pruebas de Hipótesis sobre los Supuestos
print(f"----- Pruebas de Hipótesis sobre los Supuestos para {dataset_name} -----")
# Prueba de Shapiro-Wilk para normalidad
shapiro_stat, shapiro_p = shapiro(residuals)
print(f"Shapiro-Wilk p-value ({dataset_name}): {shapiro_p:.4f}")
if shapiro_p < 0.05:
print("Conclusión: NO sigue una distribución normal.")
else:
print("Conclusión: Sigue una distribución normal.")
# Prueba de Ljung-Box para autocorrelación
ljungbox_results = acorr_ljungbox(residuals, lags=[10]) # Probando los primeros 10 retardos
ljungbox_p = ljungbox_results["lb_pvalue"].iloc[0]
print(f"Ljung-Box p-value ({dataset_name}): {ljungbox_p:.4f}")
if ljungbox_p < 0.05:
print("Conclusión: Autocorrelación significativa.")
else:
print("Conclusión: No hay autocorrelación significativa.")
# Prueba de Breusch-Pagan para heterocedasticidad
_, bp_pvalue, _, _ = het_breuschpagan(residuals, sm.add_constant(predictions))
print(f"Breusch-Pagan p-value ({dataset_name}): {bp_pvalue:.4f}")
if bp_pvalue < 0.05:
print("Conclusión: Varianza no constante.")
else:
print("Conclusión: Varianza constante.")
# Resultados de MAPE y MSE
print(f"\nResultados de Error para {dataset_name}:")
print(f"MAPE ({dataset_name}): {mape:.4f}")
print(f"MSE ({dataset_name}): {mse:.4f}\n")
# Calcular los residuales de entrenamiento, validación y prueba
resid_join_train = y_train[5:] - train_join
resid_join_val = y_val[5:] - val_join
resid_join_test = y_test[5:] - test_join
# Análisis de los residuales para cada conjunto
analyze_residuals(resid_join_train, y_train, train_join, dataset_name="Entrenamiento")
analyze_residuals(resid_join_val, y_val, val_join, dataset_name="Validación")
analyze_residuals(resid_join_test, y_test, test_join, dataset_name="Prueba")
----- Pruebas de Hipótesis sobre los Supuestos para Entrenamiento -----
Shapiro-Wilk p-value (Entrenamiento): 0.0000
Conclusión: NO sigue una distribución normal.
Ljung-Box p-value (Entrenamiento): 0.7537
Conclusión: No hay autocorrelación significativa.
Breusch-Pagan p-value (Entrenamiento): 0.0000
Conclusión: Varianza no constante.
Resultados de Error para Entrenamiento:
MAPE (Entrenamiento): 0.0252
MSE (Entrenamiento): 0.0433
----- Pruebas de Hipótesis sobre los Supuestos para Validación -----
Shapiro-Wilk p-value (Validación): 0.0000
Conclusión: NO sigue una distribución normal.
Ljung-Box p-value (Validación): 0.0151
Conclusión: Autocorrelación significativa.
Breusch-Pagan p-value (Validación): 0.0000
Conclusión: Varianza no constante.
Resultados de Error para Validación:
MAPE (Validación): 0.0226
MSE (Validación): 0.0077
----- Pruebas de Hipótesis sobre los Supuestos para Prueba -----
Shapiro-Wilk p-value (Prueba): 0.0000
Conclusión: NO sigue una distribución normal.
Ljung-Box p-value (Prueba): 0.0700
Conclusión: No hay autocorrelación significativa.
Breusch-Pagan p-value (Prueba): 0.0000
Conclusión: Varianza no constante.
Resultados de Error para Prueba:
MAPE (Prueba): 0.0359
MSE (Prueba): 0.0570